%% Specifying the conditions
% Comments and changes made by Selina in May 2020
% Get the current working directory
% based on Timmer, Marcel P., Abdul Azeez Erumban, Bart Los, Robert Stehrer, and Gaaitzen J. de Vries.
% % “Slicing Up Global Value Chains.” Journal of Economic Perspectives 28, no. 2 (May 2014): 99–118. https://doi.org/10.1257/jep.28.2.99

%data from Timmer, M. P., Dietzenbacher, E., Los, B., Stehrer, R. and de Vries, G. J. (2015), "An Illustrated User Guide to the World Input–Output Database: the Case of Global Automotive Production" ,
%  Review of International Economics., 23: 575–605, Accessed May 19, 2020, https://www.rug.nl/ggdc/valuechain/wiod/?lang=en 
workDir = pwd;

disp(['Current directory: ', workDir])
% Find the position of the string "jpe" in the workDir
jpePos = strfind(workDir, 'jpe');

% Shorten the workDir to include only up to and including "jpe"
if ~isempty(jpePos)
    workDir = workDir(1:jpePos + length('jpe') - 1);
else
    error('The string "jpe" was not found in the current working directory.');
end

disp(['Shortened directory: ', workDir])

% Define the path to the WIOD_2016 directory
wiodDir = fullfile(workDir, 'datasets', 'import', 'WIOD_2016');

q = 'shares'; 
fdc = 'ALL'; 
var = 'TOT'; 
agg = fullfile(workDir, 'datasets', 'import', 'WIOD_2016', 'agg_foreign_shares_modified_WIOD16.xlsx'); % specify where the aggregation excel file is located
pre_def='WIOD_16'; 
post_def='WIOD_16'; 
store = 'y'; % define 'y' if you want to store the final result to an excel file or 'n' if not
OUT = fullfile(workDir, 'datasets', 'import'); % specify the folder where you want to store the results
info = 'foreign shares mfg & non-mfg WIOD 2016'; % add any info between '' that you want to give to the exported excel file such that you can recognize it easily

% load data
for y = 2000:1:2014 % all years are considered (loop runs until the end of this program)

        
    c = 44; % number of countries Selina: Last country in WIOT is Rest of the World (Row) - also used
    n = 56; % number of industries
    fd = 5; % number of final demand categories
    %%Selina: Is not needed I guess - because we do not use factors here
    %load 'Foct2012'; % Employment shares and quantities, 1995,...,2009    

    %Selina: Those file directly downloaded from WIOT webpage
    clear WIOD16 

    fileName = sprintf('WIOT%d_Nov16_ROW.xlsb', y);
    filePath = fullfile(wiodDir, fileName);


    WIOD16 = xlsread(filePath, 'E7:CYK2478');
    A = WIOD16(1:n*c,1:n*c); %intermediate-input block for year y
    D = WIOD16(1:n*c,1+(n*c):(n*c)+(5*c)); % final demand block for year y
    x = WIOD16(1:n*c,1+(n*c)+(5*c)); % column vector of domestic industry gross output for year y
    va = WIOD16(6+n*c,1:n*c)'; % column (row?) vector of domestic industry value added for year y


    fd_v = D * ones(220,1); %Selina: generates the sum of the final demand across the rows (220 columns) for each country and industry
    coef = A*invd(x); % intermediate input coefficients = Selina: how much intermediates are needed to produce a unit of output of a given product ; column of intermediate-input of industry in country over total output (output at basic prices) in industry in country, "Technologiematrix: https://www.youtube.com/watch?v=YdFON0-Je6U - wieviel konsumiert �ber wieviel produziert = wieviel technolgie man braucht"   invd is a package from them - can be found in the source folder
    L = inv(eye(2464)-coef); % leontief inverse  %Selina: eye=Identity matrix - rest how one computes the leontief inverse 

    % for shares decomposition
    if strcmp (q,'shares'); %Selina: compares two strings

        if strcmp (var,'TOT');
            FINAL = diag(div(va,x))* L * diag(fd_v);  %Selina: div is a package from them - divides va by x; is in the online appendix k=F(I-B)^-1 *z
            %Selina: A typical element in k indicates the value added in the
            %production of the first final good by the particular production
            %factor (here all production factors together) employed in country
            %i and industry j
        end
        
    end



    % pre-multiplication
    if strcmp (pre_def,'WIOD_16'); % decide which country and industry groupings along the rows to load - Selina: only this used
        [pre_c] = xlsread(agg,'pre','c3:c46');
        [pre_i] = xlsread(agg,'pre','m3:m58');
    end



    c_pre_groupings = unique(pre_c(:,1)); % distilling the number of country groupings
    i_pre_groupings = unique(pre_i(:,1)); % distilling the number of industry groupings
    pre = zeros(size(c_pre_groupings,1)*size(i_pre_groupings,1),2464); % pre-specifying a matrix that has dimension country*industry by 2464

    % filling the pre-vector or pre-matrix
    for c=1:44 % all countries (includes RoW)
    for i=1:56 % all industries
    for cg = 1:size(c_pre_groupings,1) % loop across number of country groupings
    for ig = 1:size(i_pre_groupings,1) % loop across number of industry groupings
        if pre_c(c)==(c_pre_groupings(cg,1)) && pre_i(i)==(i_pre_groupings(ig,1)) % check whether country and industry grouping match
            pre(((cg-1)*(size(i_pre_groupings,1)))+ig,(c-1)*56+i)=1; % replace pre-specified matrix with a 1 if country and industry grouping match                     
        end
    end
    end
    end
    end

    % post- multiplication
    % Selina: New use this to get mfg and non-mfg industries
    if strcmp (post_def,'WIOD_16'); % decide which country and industry groupings along the columns to load
        [post_c] = xlsread(agg,'post','c3:c46');
        [post_i] = xlsread(agg,'post','m3:m58');
    end


    c_post_groupings = unique(post_c(:,1)); % distilling the number of country groupings
    i_post_groupings = unique(post_i(:,1)); % distilling the number of industry groupings
    post = zeros(2464,size(c_post_groupings,1)*size(i_post_groupings,1)); % pre-specifying a matrix that has dimension 2464 by country*industry

    % filling the post-vector or post-matrix
    for c=1:44 % all countries (includes RoW)
    for i=1:56 % all industries
    for cg = 1:size(c_post_groupings,1) % loop across number of country groupings
    for ig = 1:size(i_post_groupings,1) % loop across number of industry groupings
        if post_c(c)==(c_post_groupings(cg,1)) && post_i(i)==(i_post_groupings(ig,1)) % check whether country and industry grouping match
            post((c-1)*56+i,((cg-1)*(size(i_post_groupings,1)))+ig)=1; % replace pre-specified matrix with a 1 if country and industry grouping match
        end
    end
    end
    end
    end

    % Generating the required aggregation

    RESULT = (pre * FINAL * post); % pre and post-multiplication of the pre-specified global value chain matrix for aggregation

    test = ones(44,44);
    for i = 1:44
        test(i,i)=0;
    end
    test = kron(test,ones(2,1))'; % for getting rid of the domestic part %Selina: Kronecker tensor product between two matrixes, changed from 15 to 2 by Selina
    disp(['year is ', num2str(y)]);
    disp('Dimensions of RESULT:');
    disp(size(RESULT));
    disp('Dimensions of test:');
    disp(size(test));

    VA_for = RESULT .* test; %Selina: Element-wise multiplication 

    sumVA = sum(RESULT,1); %Selina: take the sum 
    sumVA_for = sum(VA_for, 1); 

    ratio =sumVA_for ./ sumVA;

    f_share((y-1999),:)= ratio; %%Selina: same as ratio but but in right column for the year 

    clear WIOD16

end
disp('done');

share = f_share';

if strcmp (store,'y');  % storing country results, selected industry
    filename = fullfile(OUT, [info, '.xlsx']);
    [num format] = xlsread(agg,'format','a1:q90');
    clear num;
    xlswrite(filename,format,'share','a1');
    xlswrite(filename,share,'share','c2:q90'); % save result in sheetname that can be traced back to the criteria  

end

disp('saved');
exit(0);

